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1.0 INTRODUCTION 


To provide a capability of highly accurate position and velocity 
measurements during the approach and landing phase of VTOL aircraft, NASA 
personnel at Langley Research Center are developing an RF Multilateration 
System. The system uses an angle-modulated ranging signal to provide both 
range and range rate measurements between an aircraft transponder and 
multiple ground stations. Range and range rate measurements are converted 
to coordinate measurements and the coordinate and coordinate rate information 
is transmitted via an integral data link to the aircraft. 

The objective of the work described in this report has been to conduct 
error analyses of the planned multilateration system and to determine and 
recommend data processing techniques that will provide highly accurate 
position and velocity measurements within the constraints imposed by 
computational speed and capacity. The initial studies have been concerned 
with the investigation of various proposed processing techniques and an 
investigation of errors associated with each of the techniques studied. In 
developing the processing techniques, it was necessary to continuously keep 
in mind the limitations of the ground computer. Cost and portability 
requirements limited the complexity of the computations and consequently 
also limited the computational scheme to be used. 

The following sections of the report provide a brief system description, 
describe investigations of various data processing techniques that have been 
studied, and point out advantages and disadvantages of each. For the 
techniques considered, error calculations are provided that permit a 
comparison between the various techniques. A recommended data processing 
technique is developed and error characteristics of this technique are 
discussed in detail. 

Finally, a ground-based mini-computer system is recommended which has 
both the capability of performing the recommended computations and the 
flexibility of performing additional computations (such as Kalman filtering) 
if future requirements make such computations desirable. 
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2.0 SYSTEM DESCRIPTION 


The planned multilateration system configuration is as shown 
generally in Figure 2-1. A single transmitter transmits a ranging signal 
to an aircraft transponder where the signal is offset in frequency and 
retransmitted to four or more ground receivers. The ground receivers 
contain digital processors which provide digital outputs proportional 
to the slant range between the receiver and the transponder and the 
rate of change of slant range. The digital receiver outputs then go to 
a central ground data processor which transforms the slant range measure- 
ments to X, Y, Z position and rate data. In the operational system, the 
position and rate data are then transmitted to the aircraft over a data 
link at the same rf frequency as the ranging signal and are available at 
the aircraft for use in navigation and guidance equipment. 

The location of the ground receivers is chosen to provide high 
accuracy throughout the touchdown area. For situations in which it is 
desirable to have accuracy over a 360° azimuth range, a configuration such 
as that shown in Figure 2-2A appears desirable. The central receiver is at the 
touchdown point. For situations in which there is a preferred direction 
of approach, a configuration such as the one shown in Figure 2-2B may be more 
desirable. To provide higher accuracy over certain segments of the approach 
trajectory, additional receivers may be used. 

The functional operation of the system may be understood by referring 
to Figure 2-3. The transmitter generates a ranging carrier, a phase modulated 
tone for ranging, a reference carrier, and a data subcarrier. At the 
aircraft transponder, the incoming signals are translated in frequency 
coherently using the modulation on the reference carrier. The ranging 
carrier and associated tone ranging modulations are then retransmitted to 
the ground receivers. At the ground receivers, the Doppler shifted carrier 
and ranging modulation are processed to derive digital signals proportional 
to the two-way time delay (range) and the rate of change of time delay 
(range-rate) . These digital measurements are then transmitted to the coordinate 
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(B) Approach from North only 


Figure 2-2. Plan view of possible ground receiver locations for four 
ground receivers. 
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Figure 2-3. Simplified functional block diagram. 







computer for derivation of coordinate position and coordinate rates. 

The theoretical operation of the system may be seen by examining the 
mathematical expressions for the signals throughout the system. The trans- 
mitted signal may be written as 

v fc = V t ^sin[N^a) Q t + gsinat] + V^sint^a^t + ysin w Q t] (2-1) 

where 

= multiplication factor (1760) 

N£ = multiplication factor (1764) 
aj Q = reference oscillator frequency (3.1 MHz) 
a = ranging tone frequency (100 Hz) 

6 = ranging tone phase deviation ( 5280 rad.) 
y « reference carrier phase deviation . 

At the transponder in the aircraft, the signal received is similar to 
equation (2-1) but is delayed by a time delay t^. At the transponder, 
the signal is shifted in frequency using the reference signal, and is 
retransmitted and received at a ground receiver after an additional 
delay t 2 as 

v R = V R sin[N x m o (t - - t 2 ) + 3sin a(t - - t 2 )] (2-2) 

where 

= uplink delay 

t 2 = downlink delay 

N = N, - 32 
x 1 

and the other nomenclature is defined above. 
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At the receiver > the incoming signal is mixed with a local reference: 


V REF = V REF ® ln[ »i“o t + &Sin ° t] 


(2-3) 


and the resulting first IF signal is given by 


V IF1 = V IF1 [(N 1 ‘ V V + N x“o (T l + t 2 } ‘ 6(sin at 


- >sin a(t - (x^ + x 2 )) ] . 


(2-4) 


An additional mixing operation takes place with a multiplied version of 
the reference oscillator signal to provide a second IF signal given by 


V IF2 ■ v if2 sln[ V + \"o <T ! + T 2> 


a ^ T l + x 2^ (t. + T_) 

+ 23(sin = cos a (t - 


)} (2-5) 


where the form of the last term has been modified using a trigonometric 
identity for the difference between two sinusoidal signals. 

The instantaneous frequency of the IF signal may be determined 
by differentiating the argument of the sine function above and is 
given by 

(t, + X 2^ 

03 =03 + N 03 (x n + t„) + 23a[sin a(t r ) 

IF o x o 1 l l 


sin 


°^ T l + 


] + $ax cos a(t - x n - x 2 ) 


(2-6) 


where t is the rate of change of the propagation delay. The last term in 
equation (2-6) is an error term and system parameters are chosen to make this 
term small in comparison to the desired terms in the expression. Thus, the 
IF frequency is given approximately by 




Figure 2-4. Sketch of IF frequency showing range and range rate (Doppler) terms. 


to 


IF 


to + N to (i. + t„) + 26a sin 
o x o 1 2 


/ 0(1^ + t 2 ) 
\ 5 


sin a(t - 



(2-7) 


A sketch of this equation is shown in Figure 2-4. As shown on the sketch, 
the deviation of the IF frequency over a modulation cycle is proportional to 
the total time delay or two-way range, and the average IF frequency is equal 
to the reference oscillator frequency plus a Doppler offset proportional to 
the two-way rate of change of delay, or range rate. 

A technique of cycle counting is used to obtain measures of two-way 
range and range rate. Three gated counters are used to count the IF 
frequency over the first and second halves of the modulation cycle and the 
reference frequency over a full modulation cycle. The sum of the counts 
over the first and second halves of the modulation cycle is proportional to 
the average frequency and the difference in these counts is proportional to 
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the frequency deviation. The reference frequency counter is used to obtain 
the difference between the average received frequency and the transmitted 
frequency. 

The number of counts over the first half cycle is given by 


7r/a + t/2 



t/2 


where e^ is a fractional number 0 _< e^ < 1 selected to make the number of 
counts* C, be an integer value. In the integration, t = is assumed 

constant over the duration of the modulation cycle. Over the second half 
cycle. 


2ir/a + t/2 



ir/a + t/2 


where 0 < < 1. The reference counter gives 


(2-9) 


2ir/a 

C_ = ~ /a) dt - e. 




( 2 - 10 ) 


where 0 < e^ < 1. 

Evaluation of the integrals using equation (2-7) and summing and differ- 
encing the counters gives, for the range rate count, 


C R " C 1 + C 2 C r “ a (x l + X 2 ) -e l e 2 + e 3 


( 2 - 11 ) 


md for the range count 


C = C - C — — - 
C R C 1 2 it 


t + T 

sino( = ) 


e l + e 2 


( 2 - 12 ) 




where (-2 < < 1) and (-1 < e R < 1) . 

Thus, the counters provide digital outputs proportional to the two-way 
range and range rate. These digital outputs are then used to calculate the 
one-way ranges and range rates and the coordinate position and velocity 
components of the aircraft. Conversion of the range and range rate measure- 
ments to coordinate position and velocity is accomplished in the ground 
computer. 
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3.0 SOLUTION TECHNIQUES 

3.1 General Discussion 

As seen in the previous section, the basic system measurements are the 
range and rate between the ground stations and the airborne transponder. 

These basic measurements will contain random and bias errors due to 
system characteristics, quantization errors, and calibration errors. It 
has not been within the scope of the present work to investigate the 
sources of these errors in detail, but to consider only the effect of the 
basic measurement errors on the solution for coordinate velocities and 
position. The central data processor in the system has the function of 
converting the basic measurements to coordinate rates (x, y, z) and 
coordinate position (x, y, z) . These coordinate positions and rates are 
then transmitted to the aircraft for use in the navigation and guidance 
system. The basic measurements will be taken at a data rate of approximately 
ten per second. 

The overall accuracy achievable with the multilateration system will 
depend strongly on the data processing techniques used. For example, with 
four or more ground stations, a redundant measurement is available which 
may be combined with other measurements to improve the system accuracy. 

It is also possible to provide in-flight calibration so that a sequence of 
measurements is used to correct for bias in the system during the initial 
stages of an approach path. In general, it is desirable to use all of the 
basic measurement information available to provide the most precise position 
and velocity estimation possible. 

Since a fast-time response is desired in the system, it is preferable 
if possible to use point estimation techniques; that is, to calculate 
coordinate position and velocities each sample instant based on the 
measured data. While the accuracy can be Substantially improved by using 
the time history of the measurements in a filtering technique (e.g., 

Kalman filtering), these techniques have not been considered in this study. 
They will, however, be considered in future work. 
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Fig. 3-1. Definition of vectors. 



The equations necessary for conversion of range measurements to 
coordinate positions are non-linear, hence the determination of the best 
solution technique is not a straightforward problem. In this report, 
several solution techniques have been developed and the advantages and 
disadvantages of each technique have been determined. For convenience, 
the solution techniques developed have been broken down as follows: 

A. Solutions for position coordinates. 

1. Explicit solutions for three ground stations. 

2. Linearized minimum variance solution. 

3. Iterative minimum variance solution. 

4. Explicit minimum variance solution for N ground stations. 

B. Solution for coordinate rates. 

1. Explicit minimum variance solution for N ground stations. 

C. Bias (system time delay) removal techniques. 

Solutions combining rate estimation with position estimation have not 
been considered because of the computational complexity involved. 

The following sections discuss the above solution techniques and outline 
the considerations that lead to the determination of the recommended solution 
technique discussed in Section 6.0. 

3.2 Explicit Solutions for Aircraft Position 
Coordinates with Three Ground Stations 

In vector terminology, the range vector from a receiver i to the 
aircraft transponder may be written as 

R. = R - (3-1) 

where the vectors are defined in Figure 3-1. The magnitude squared of the 
vector iL is 


2 2 
R* » R + 


P? - 2R 


-2 


(3-2) 
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Thus, one equation of the above form will be obtained for each ground 
receiver in the system. A minimum of three stations is required to solve 
for the three unknown position coordinates of the aircraft (x ,y ,z ) where 
is assumed positive in all cases and the station locations are known. 

For three stations, the coordinate system can always be chosen so 
that one station defines the origin of coordinates, one station defines the 
X axis., and the third station defines the XY plane. Therefore, let the 
station locations be as follows. 


Station 

1 

2 

3 


Coordinates 
( 0 , 0 , 0 ) 
(0, y 2 , 0) 
(x 3 , y 3 , 0) 


With this choice of coordinate system, the following equations apply: 


R I = 

R 2 * 

R? = 


2 . 2 
x + y + z 
a a 


2 

x + y + z 
a J a a 


2 2 2 
a + *3 + y 2 ' 2 V2 

2 


2 2 2 2 z . 

x +y +z + x_ + y„ - 2 x x_ 
a ■'a a 3 ■'S a 3 


- 2 y a y 


a' 3 


(3-3) 

(3-4) 

(3-5) 


and the solution for x , y , z is straightforward. First solve for y , 

3 3 3 3 


y, ■ 


a 2y, 


(r; - 


n 2 A 2. 
r 2 + y 2 ) 


(3-6) 


then solve for. x , 
a 


1 2 2 2 2 
x a = 2x^ (R 1 - r 3 + x 3 + y 3 


2 V 3 > 


(3-7) 


and finally, solve for z r 


/D 2 2 2 . 1/2 
z = + (R, - x - y) 
a 1 a 7 a 


(3-8) 
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This solution in the coordinates x, y, z can then be obtained in 
another coordinate system X, Y, Z by using a linear transformation 

X = [K] x (3 “ 9) 

where 

x = 

and [K] is a 3 x 3 matrix for the linear coordinate transformation. 

A completely general explicit solution for the three station case in 
an arbitrary coordinate system is derived in Appendix A. A technique for 
error analysis of this type of solution is given in Appendix B. 



3.3 Linearized Minimum Variance Solution for N Ground Stations 

3.3.1 The weighted least squares technique .— With more than three 
ground receivers, redundant information is available which may be used to 
improve the accuracy of the aircraft position estimate. For example, with 
four ground receivers, four estimates of aircraft position may be obtained 
using all combinations of three of the equations relating the measurements 
to the position coordinates. To combine the measurements in some optimal 
manner, the technique of weighted least squares estimation is used. This 
technique provides a means of linearly combining solutions so that the 
resultant estimate has minimum variance. 

To demonstrate this technique, consider for example two measurements, 
y^^ and y 2> of a parameter y and of uncertainty implied by their 
respective known standard deviations and To determine an estimate 

of the parameter y given the measurements, we form a linear weighted average, 

y = W 1 y l + w 2 y 2 (3-10) 

where w^ and w 2 are weighting factors and w^ + w 2 = 1. 

A 

The variance of errors associated with the estimate y is calculated from 
the above equation as 
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(3-11) 


2 2 2 2 2 
°y = W 1 °1 + W 2 °2 

where, for simplicity, y 1 and y_ are assumed independent with zero 

l / 2 

correlation. To determine the w. that minimizes the variance a-', dif- 

i y* 

ferentiate eq. (3-11) with respect to w x and set the result equal to zero. 

2 2 2 2 2 

°y = W 1 °1 + ^ " w i^ a 2 (3-12) 


^^°y ^2 2 

' J W Z 1 ' = 2 w i a i + 2(1 - w x )(-l) o 2 = 0 . 


(3-13) 


Therefore the weighting factors are 


w, = 


2 2 
°1 + °2 


"2 „ 2 2 
°1 + °2 


and 


2 2 
a l + °2 


min variance 


y i + 


2 2 
°1 +a 2 


(3-14) 


The variance of y is minimized and is given by 


o* 

y 


2 2 
°1 a 2 
2 2 
0 1 +0 2 


min 


(3-15) 


The above result can also be obtained by selecting y to minimize the 
weighted mean square error risk function R, where 
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(3-16) 


R = 



Minimizing with respect to y gives 



y) 



y) 



0 



square error 


(3-17) 


(3-18) 


The matrix formulation of the weighted least squares or minimum variance 
technique provides results as shown above in a very compact notation. In 
ref. 1 it is shown that if the measurement vector y is related to the 
estimator vector b by an equation 

y = [A] b + v (3-19) 

where v is the (zero mean) measurement noise vector, then the estimate vector 
which minimizes a weighted least squares risk function is 

b = (A T tp' 1 A) -1 A T 4> _1 y (3-20) 

where is the weighting matrix. For minimum covariance, the weighting 
function ip must be the covariance matrix of the measurement errors (^-) . 

If this is the case, the covariance matrix of the estimate error is 

- < aT a)_1 ' <3 ' 21) 

The matrix technique can be demonstrated by applying it to the example 
previously discussed. For the two measurements y^ and y 2 we have 
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y l ■ y + V 1 
y 2 - y + v 2 

Cov 

Thus in this case 



2 

Var Ay^ = Var v^ * a ^ 

2 

Var Ay 2 = Var v 2 - 
Ay^, Ay 2 = Cov v^v 2 - 0 . 



(3-22) 

(3-23) 

(3-24) 


and the answer follows directly from eq. (3-20). The various calculations in 
obtaining the solution are as follows: 



(3-25) 


(3-26) 




(3-29) 


which is the same equation obtained previously. 

The advantages of the matrix notation in dealing with four or more 
equations and correlated measurement errors should be evident. 
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3.3.2 Linearized solution for the special case of all ground stations 
in a plane . — For all ground receivers in a plane, the range equations (see 
Figure 3-1 for definition of vectors). 


- 2R • i = 1, n t (3-30) 

2 2 
can be linearized by treating R as an independent variable. Letting u = R /2 

and 

_2 2 
K - Pi 

q ± = 2 = “ x i x - y ± y + u (3-31) 

gives the equation 


Q = [Pi X (3-32) 

where, for n =* 4, 



The solution is given by 

x = (P T v"J P)' 1 P T ^ Q (3-33) 


where ’{' is the covariance matrix of errors in the q.. Since 
AQ i 


Aq ± z R i AR ± and AQ = [T] AR 


(3-34) 
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we have 


% • E l"! *5*1 - 1 tT 


(3-35) 


where is the covariance matrix of the slant range measurements (derived in 
Appendix C) and 

^000 


T = 


0 r 2 0 0 

0 0 r 3 0 


0 R 


4 -* 


After calculation of x, y, and u the z estimate is found as 


*2 * 2 . 1/2 
(2u - x - y ) 


(3-36) 


An error analysis of this solution technique is given in Appendix D. 

3.3.3 Linearized solution for the special case of all ground stations 
not in a plane . — For four or more ground stations not in a plane, a linear 
solution can be found using range difference equations. Using the same 
nomenclature as before, define 


2 _2 , 2 2 
R i - \ + P n - Pi 


(3-37) 


By subtracting the nth range equation from the remaining equations, we 


obtain 


= - R * (P ± - P n ) (3-38) 

or in matrix form 

Q = [P]* . (3-39) 
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where 



[P] 


x,-x 
1 n 


y r y n 


z,-z 
1 n 


x 0 -x 
2 n 


y 2 -y n 


z 0 -z 
2 n 


s i -x y , -y 
n-1 n n-1 7 n 


z ,-z , 
n-1 n 


x = 


x 

y 

z 


For four ground stations not in a plane, the solution is unique and is 
given by 


x = (p T p) 1 p T Q . (3-40) 

T -1 

Note that if the stations lie in a plane, the inverse (P P) does not exist. 

For more than four ground stations not in a plane, the weighted least 
squares solution may be used. This solution is 


x 




Q 


(3-41) 


where is the covariance matrix of the errors in the vector Q (see Sec. 3.3.2). 
An error analysis of this solution is given in Appendix E. 

3.4 Iterative Minimum Variance Solution 


The derivation of an iterative minimum variance solution is given in 
Appendix F. The technique is summarized as follows. 

The two-way range measurements (R^ + R^) are related to the coordinate 
positions as follows: 


A ± = [(x - x x ) 2 + (y - yi ) 2 + (z - zp 2 ] 172 + [(x - x ± ) 2 + (y - y.) 2 


. , . 2 . 1/2 

+ (z - z ± ) ] 


(3-42) 


i = 1, n , 
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or 


A i = f i^ x * y » • 


The above equations are expanded in a Taylor's series around an initial 
estimate of the solution to obtain 


i 9£ n 9£ n 9£ il 

£ 1 U,y, Z )J o + Sr\ 0 Lx+ 3rJ„ ay+ 3rJ n 4 ° ' A i ' (3 ‘ 43> 


Thus , 


Mt-Ai 


1 5£ il 9£ il 9£ ll 

- Ax + 35-J o Ay + -jj-J 


Az 


(3-44) 


The set of equations above can be written in matrix form as 


AA = [D] Ax 


(3-45) 


where 


AAl 

AA = ) 


AA 


Ax = 


Ax 

Ay 

Az 


and [D] is a n x 3 matrix of partial derivatives of the f^. 

Solving for Ax from eq. (3-45) gives the least squares form, 


T -1 T 

Ax = (D D) D AA 


(3-46) 


and an improved estimate of x is obtained over the initial estimate where 

\+i ■ \ + • (3 ' 47) 

This improved estimate is then used as the initial solution and the process 
is repeated. 

Thus, the iterative solution procedure is : 

1. Determine initial solution by some means. 

2. Calculate D matrix. 
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3. Calculate AA. 

4. Calculate Ax. 

5. Calculate improved estimate of x. 

6. Go to step 2. and repeat calculations until | Ax | is smaller than 
some predetermined value. This indicates the solution has 
converged to the value for x within the limit preselected. 

The covariance matrix for the estimation errors is, from eq. (3-46), 


ip__ = E{Ax Ax T } = (D T D) 1 D T <p_ D(D T D) _1 . (3-48) 

Ax AA 


3.5 An Explicit Minimum Variance Solution for n Ground 
Stations in an Arbitrary Coordinate System 

For n ground stations, n equations of the form 

R ± 2 = R 2 + Pi 2 - 2R • p t i - 1, 2 n (3-49) 


are obtained where R^ is the measured slant range and p^ is the vector 
from the origin of coordinates to station i. R is the slant range vector 
from the origin to the aircraft. 

Now, subtracting the n^ equation from the other n-1 equations gives 


i 

n r i 

2 

rn = R • 

(5 n - p 4 > 

i— 1, 2, . . • , (1 n) (3—50) 

or 


q t = R * 

<p n - PP 

(3-51) 


where q^ equals the term on the left-hand side of equation (3-50) . 

The set of equations (3—51) can be expressed in matrix form by 

Q = [P]x (3-52) 
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where 


Q - 


4 2 
q n-l 


x = 


x 

y 

z 


p = - 


X 1 

- X 

n 

y i 

' y n 

Z 1 


z 

n 

X 2 

- X 

n 

y 2 

- y n 

Z 2 


z 

n 


• 


• 


• 



• 


• 


• 



• 


• 


• 


‘n-1 

- X 

n 

y n-l 

- y n 

Z n-1 


z 

n 


(n-1) x 3 


Now partition [P] and x as follows: 
Q- IAB]{’} 


where 


B = - 


z. - z 
1 n 


L Z n-l ' Z nJ 


(n-1) x 1 


and 


A - - 


X 1 " x n 
x 2 " x n 


X - - X 

n-1 n 


y l - y n 
y 2 " y n 


y n-l " y n 


(n-1) x 2 


(3-53) 


With this partitioning, we have 

Q = Ap + Bz 

To obtain a least squares estimate of p, we assume that Q - Bz is the 
observed data. Under this assumption. 


(3-54) 


p *.s. = CQ " CBZ 


(3-55) 
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where 


(3-56) 


C - (A%'J A)' 1 . 


The expression is the covariance matrix of errors in (Q - Bz) . 


(Note: For all ground stations in a plane, B = 0 and equation (3-44) is 

considerably simplified, 
the covariance of errors in Q only.) 

The original equations (3-49) may be written in terms of p as follows: 


In this case, the covariance matrix ip becomes 

gB 


Rj. 2 = (P - Pj) 1 (P ~ P ± ) + z 2 - 2z^z + z 2 i = 1, 2, 


n 


(3-57) 


where 



A 

Any of the n equations may be used to solve for z by substitution of 
the estimate for p, equation (3-55), into equation (3=57). When this is 
done, the resulting equation can be solved for z to yield: 


z(i) = 


-b . +■%/ b. 2 - 4ac . 

i — J l i_ 

2a 


i = 1, n 


(3-58) 


where a 



c. 

i 


1 + b t c t cb 

-2z i - 2Q T C T CB + 2p i T CB 

- R i 2 + p ± T p i - 2p J T CQ + Q T C T CQ + z A 2 . 
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The z estimate z can be obtained as a weighted average of the n 
z(i) values obtained above. Thus, we have for 2 


z = (U 1 ^ 1 U)- 1 U^ 1 z 


(3-59) 


where 


U = 


1 


z(l) \ 

1 1 


z(2) / 

• g 

1 

j 1 
■ 

z ■ ! 

1 • 

n x 1 | 

i z(n) ) 


and ip is the covariance matrix of the errors in z(i). This matrix is 
derived in Appendix G. The standard deviation of the z estimation error is 


2 T -1 -1 

°A2 • = (u \z U > 


(3-60) 


The x and y estimates now determined from eq. (3-55) are 




CQ - CBz . 


(3-61) 


Error covariances associated with this solution technique are derived in 
Appendix G. 


3.6 Minimum Variance Solution for Coordinate Rates 
in an Arbitrary Coordinate System 

For the velocity components, the measurements made by the receivers 
are (for four receivers) 

M i « R x + R ± i = 1, 2, 3, 4 (3-62) 

where R is the range rate between the transponder and receiver i. In 
vector terminology, the magnitude of the range rate is given by 
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(3-63) 


X 



where x = 1 x+1 y + 1 z and R. is the range vector from receiver i to 

x y J z i 

the transponder. Since the estimate of can be calculated first, the 
velocity estimate can be obtained from 


a ♦ 

R • x 


A i = 


i = 1, 2, 3, 4 . 


where 


• • 


A ± = M i - M^^/2 


i = 1, 2, 3, 4 . 


In matrix form. 


where 


A * [S] x 


A = 


4x1 


x 

y 

z 


3x1 


(3-64) 


(3-65) 
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and the solution is given by 


- T -1 -1 T -1 - 

* - < s *> s *sA A 


( 3 - 66 ) 


where is the covariance matrix of the range rate measurement errors. 

The explicit form of this matrix is derived in Appendix C and, for the case of 
equal error variance in each two-way measurement, is given by (for four 
ground stations): 



where a,* is the variance of the two-way range rate measurement error. The 
AW 

error covariance matrix for coordinate rates is derived in Appendix H. 
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3.7 Time Delay Bias Removal Techniques 


3.7.1 General technique .— The measured ranges from the ground receivers 
to the aircraft transponders will contain bias due to unknown time delays in 
the transmitter, transponder, and receivers. For the four-station case, 
measurements (M^) made by the four receivers are 


= R 1 + R ± + i = 1, 4 (3-68) 

where is an error corresponding to an uncompensated delay in the R^, 

R^ path. 

A suggested technique for removal of the uncompensated bias errors 
involves the estimation of the error from the equation 

e i = M i ' “i » (3-69) 

where = R^ + R^ is the calculated estimate of the two-way slant range. The 
estimate is then smoothed by using a long time constant low-pass digital 
filter such as 


„ k+1 




(3-70) 


where t is the filter time constant and At is the calculation increment. 

A flow chart of this calculation technique is shown in Figure 3-2 . 
3.7.2 Variance of errors in time delay bias estimation .— The variance 
of errors associated with this estimation technique can be determined as 
follows. A small change in the estimate is written 


Ae ± = AM i - 



(3-71) 
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/s ^ 2 2 2 1/2 * 

Also, a small error in = [ (x - x^) + (y - y^) + (z - z^) ] + may be 

written In terms of errors in x, y, z 


9M ^ 9M. 9M 

AM. = — ^ Ax + — Ay + — - Az (3-72) 

1 9x 9y 9z 


or in matrix form 


AM, = [Ej Ax 
i i 


where 


and 


Thus , 


E. 

l 


3M 3M. 3M. 

i x l 

A 5 > 

3x 3y 3z 


Ax 



As. = AM. - E,Ax . 
1 X X 


(3-73) 


(3-74) 


An upper bound on the bias estimation error can be determined by 
assuming unity correlation between errors in an individual measurement M 
and errors in the estimated value of the measurement M. This upper bound 
is given by 


Ae i < 


+ { e i e i t } 


(3-75) 


where is the variance of the two-way range measurement and is the 

covariance matrix of the estimation error. The exnression ilrr - is derived in 

T Ax 

the appendices for the various solutions considered. 
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A more exact estimate of the bias error estimation can be obtained 
if the relationship between the two-way measurement errors AM and the Ax 
is known. This relationship depends on the solution technique used; for 
example, for the iterative method described in Appendix F, we have 

Ax = G -1 D T AM (3-76) 

where AM in the above equation corresponds to AA in the nomenclature used 
in the appendix. Using eqs. (3-74) and (3-76) we have 

Ae ± = AM t - E i G _1 D T AM (3-77) 


or in another form 

Le ± = AM - E ± G _1 D T AM (3-78) 

where 

= (1, 0, 0, 0} , 8 2 * (0, 1, 0, 0} , etc. 

Thus, the variance of a single measurement of a time delay bias is 
given by 

a L t - [ »i - e i G_1 ” T) [e i - h G ' x dT]T < 3 - 79> 

where drr— is the measurement error covariance matrix. This variance can be 
AM 

considerably reduced by filtering as is shown in the following section. 

3.7.3 Reduction of time delay bias estimation errors by filtering .— 

Using the digital low-pass filter given by eq. (3-70), the variance of errors 

- k + 1 

in the smoothed estimate is given by 

Var + 1 = (1-a) 2 Var + a 2 Var , (3-80) 
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where a = At/x. For times t >> x the smoothed variances will not change 
with time so that for t << x 

Var + 1 « Var (3-81) 

and the variance of the filtered output may be written in terms of the 
variance of the unsmoothed estimate of f 


Var e. = -r- - — Var e. ; t >> x 

i 2 - a i 


(3-82) 


A At 

or approximately, Var as Var ; t >> x >> At. Thus, for a calcu- 
lation increment of .1 second and a filter constant x of 2 seconds, the 
variance of the smoothed estimate is 2.5% of the unsmoothed estimate. Longer 
filter time constants provide correspondingly better variance reductions. 

It should be noted that instead of a conventional low-pass filter, a 
minimum variance (i.e., optimal) filter could be used. For this type of 
filter, the parameter a in the filter equation is recalculated for each time 
increment, starting with an initial value of unity. The recursive relation- 
ship for minimum variance is derived as 


k + 1 
a 


k 

a + 1 



(3-83) 


In using this technique, the value of a will eventually converge to zero, 
corresponding to an infinite time constant filter. For this reason, it is 
desirable to limit the minimum value of a to a small number on the order 
of .01 (i.e., the number .01 corresponds to a low-pass filter of 10 second 
time constant when using a calculation increment of .1 second). The advantage 
of the optimal filter is that initial values of the estimate are more 
accurate than those from a conventional low-pass filter with its inherent 
transient buildup to the steady state. 
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4.0 ERROR STUDIES 


4.1 General Discussion 

This section considers the errors associated with the various solution 
techniques discussed in the previous section. The error variances associated 
with the solution for position coordinates as a function of the slant range 
measurements are a function of the following system parameters: 

1. Basic slant range measurement errors, 

2. Station location, 

3. Aircraft location, 

4. Solution technique, 

5. Number of ground stations. 

The errors in estimation of the coordinate rates are functions of all 
the above parameters plus the additional parameters : 

6. Aircraft velocity. 

_ 7. Range rate measurement errors . . _ _ 

In order to present the results of the error studies in a useful form, 
several techniques have been used. For the position errors, contour plots 
are useful for visualizing the rms errors that will be found at various 
geometrical locations. For presenting results on errors associated with 
coordinate rates, specific trajectories must be defined and used in the 
error calculations. In this section, several specific trajectories have 
been defined and the position and rate errors associated with these 
trajectories have been calculated. These results are shown in the following 
sections . 


4.2 Comparison of Position Error Variances 
for Various Solution Techniques 

To compare the efficiency of the solution techniques discussed in 
Section 3.0, tables have been prepared to compare the solution techniques 
for specific station locations and aircraft locations. The derivation of 
the error equations is given in the appendices. 

In addition to the theoretical error variances, computer simulation 
techniques have been used to determine actual variances that might be 
encountered under experimental conditions. For these computer simulations. 
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the technique shown on the flow chart of Fig. 4-1 is used. A random number 
generator is used to calculate random perturbations on the slant range 
measurements and the various solution techniques are used to solve for the 
coordinate positions. These coordinate positions are then compared to the 
true coordinate positions to determine coordinate errors. Five hundred 
samples of this calculation are made, and the mean and variance of the 
coordinate errors are determined. This simulation technique gives confidence 
in the theoretical results in cases where the numerical values compare 
favorably. 

Table 4-1 provides a comparison of various solution techniques for an 
aircraft at the position (1000 ft, 1000 ft, 100 ft) and for the station 
locations given in the table. In this case, the position of the aircraft 
is somewhat out of the operational range of the system (i.e., below 15° 
elevation); however, an extreme case was desired to provide large errors 
for comparative purposes. As may be seen, the iterative solution, explicit 
solution and the linearized minimum variance solution provide the same 
standard deviations of coordinate errors. The three-station trilateration 
solution is the worst of the solution techniques considered, as would be 
expected. Values from the computer simulation compare favorably to the 
theoretical values, and it should be noted that the iterative solution con- 
verges quite rapidly (i.e., only a few iterations are required). 

Table 4-2 provides the same type of comparison, except that two of the 
ground stations are elevated by 100 ft. Notice that the standard deviations 
in the z coordinate are considerably improved over the values with all ground 
stations at z = 0. Standard deviations of x and y are also improved by 
elevating the ground stations. Tables 4-3 and 4-4 provide the same type of 
comparison, except for an aircraft position of (1000, 1000, 500). Notice 
that at reasonable altitudes, the standard deviation in the z coordinate 
improves dramatically. Tables 4-5 and 4-6 provide the same comparison for 
the aircraft located at (100, 100, 50). As the aircraft approaches the 
landing position, the benefit of elevated stations decreases. In all cases 
with elevated stations, the iterative solution provides the best accuracy, 
followed closely by the explicit minimum variance solution. For non-elevated 
stations, all solution techniques give the same accuracy. 

All of the above calculations were made with a fixed-station configuration 
as given in the tables. The effect of changing station locations is considered 
in the following section. 
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REPEAT FOR 
NUMBER OF 
SAMPLES 
DESIRED 


Figure 4-1 • Flow Chart of Technique for Direct Computation of Mean and 

Standard Deviation of Coordinate Errors by Computer Simulation. 
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Table 4-1. Comparison of error standard deviations for various solution 
techniques for an aircraft at the position (1000,1000,100) 
and for the station locations indicated. 


Solution Techniques and 
Station Locations 

Standard Deviations of Error Assuming 
2 ft rms Two-Way Range Measurement Error 

Station Locations 

01 ( 0, 500, 0) 

02 ("433, -250, 0) 

03 ( 433, -250, 0) 

//4 ( 0, 0, 0) 

Theoretical 

Values 

Values From 
Computer 
Simulation 

0 

X 

a 

a 

z 

0 

X 

a 

y 

0 

Z 

1. Three station (trilatera- 
tion) solution (01, 02, 03) 


4.84 

70.7 

5.18 

4.31 

45.9 

2. Ti'np.ori v a A m'mitnl 

variance solution 

5 -42 

■ 

on o 

t Vr • 

5.07 

/. on 

-T « «.V 

q q /. 

s -t • —r 

3. Explicit minimum variance 
solution 

5.42 

4.78 

70.3 

5.07 

4.20 

43.7 

4. Iterative solution 

5.42 

4.78 

70.3 

4.06 

3.80 

41.79 
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Table 4-2 . Comparison of error standard deviations for various solution 
techniques for an aircraft at the position (1000,1000,100) 
and for the station locations indicated. (2 Stations Elevated) 


Solution Techniques and 
Station Locations 

Standard Deviations of Error Assuming 
2 ft rms Two-Way Range Measurement Error 

Station Locations (ft) 
n ( 0, 500, 0) 

n (-433, -250, 100) 
//3 ( 433, -250, 100) 
#4 (0, 0, 0) 

Theoretical 

Values 

(ft) 

Values From 
Computer 
Simulation 
(ft) 

a 

X 

H 

a 

z 

a 

X 

a 

y 

o 

z 

1. Three station (trilatera- 
tion) solution (#1, //2 , // 3) 

5.45 

4.22 

A2A 

5.17 

4.40 

45.3 

O T i-. — -T 1 J - J w.. . 

4. • ljXUCdi.i.4<CU IUX.liJ.mUUl 

variance solution 

5.45 

7.61 

49.7 

5.08 

7.52 

49.7 

3. Explicit minimum variance 
solution 

5.83 

3.93 

36.6 

5.45 

3.67 

36.2 

4. Iterative solution 

4.69 

3.88 

34.0 

4.62 

3.61 1 

32.6 
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Table 4-3. 


Comparison of error standard deviations for various solution 
techniques for an aircraft at the position (1000,1000,500) 
and for the station locations indicated. 


Solution Techniques and 
Station Locations 

Standard Deviations of Error Assuming 
2 ft rms Two-Way Range Measurement Error 

Station Locations 
#1 ( 0, 500, 0) 
If 2 (-433, -250, 0) 

If 3 ( 433, -250, 0) 

H (o, o, 0) 

Theoretical 

Values 

Values From 
Computer 
Simulation 

0 

X 

0 

y 

o 

z 

0 

X 

a 

y 

o 

z 

1. Three station (trilatera- 
tion) solution (//l, #2 , If 3) 

5.69 

5.10 

14.7 

5.53 


15.4 

2. Linearized minimum 
variance solution 

cr r. c 

J .DJ 


B 

fl 

4.41 

15.3 

3. Explicit minimum variance 
solution 

5.65 

5.06 

14.7 

5.71 

4.41 

15.4 

4. Iterative solution 

5.65 

5.06 

14.7 

5.71 

4.41 

15.3 
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Table 4-4. Comparison of error standard deviations for various solution 
techniques for an aircraft at the position (1000,1000,500) 
and for the station locations indicated. (2 Stations Elevated) 


Solution Techniques and Standard Deviations of Error Assuming 

Station Locations 2 ft rms Two-Way Range Measurement Error 


Station Locations 

Theoretical 

Values 

Values From 
Computer 
Simulation 

01 ( o, 500 , o) 
in (-433, -250, 100) 
//3 (+433, -250, 100) 
#4 (0, 0, 0) 

o 

X 

o 

y 

O 

z 

0 

X 

a 

y 

n 

1. Three station (trilatera- 
tion) solution ( ill , 02, 03) 

5.60 

4.07 

12.9 

5.45 

3.44 

D 

2. Linearized minimum 
variance solution 

5.60 

8,11 

52,4 

5,64 

7.26 1 

50.6 

3. Explicit minimum variance 
solution 

5.69 

4.20 

12.8 

5.74 

3.49 

13.25 

4. Iterative solution 

5.57 

4.07 

12.7 

5.20 

3.50 

12.17 
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Table 4-5. Comparison of error standard deviations for various solution 
techniques for an aircraft at the position (100,100,50) 
and for the station locations indicated. 


Solution Techniques and 
Station Locations 

Standard Deviations of Error Assuming 
2 ft rms Two-Way Range Measurement Error 

Station Locations 

01 ( 0 , 500 , 0 ) 

02 (-433, -250, 0) 

03 (+433, -250, 0) 
#4 ( 0, 0, 0) 

Theoretical 

Values 

Values From 
Computer 
Simulation 

°x 

0 

y 

a 

z 

n 

a 

y 

a 

z 

1. Three station (trilatera- 
tion) solution (01, 02, Ok) 

2.48 

1.28 

5.3 


1.77 

7.34 

2 . Linearized minimum 
variance solution 

mj 

fl 

C O 

B 


7.46 

3. Explicit minimum variance 
solution 

1.81 

1.27 

5.3 

1.91 

1.69 

7.12 

4. Iterative solution 

1.81 

1.27 

5.2 

1.91 

1.69 

7.46 
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Table 4-6. Comparison of error standard deviations for various solution 
techniques for an aircraft at the position (100,100,50) 
and for the station locations indicated. (2 Stations Elevated) 


Solution Techniques and 
Station Locations 

Standard Deviations of Error Assuming 
2 ft rms Two-Way Range Measurement Error 

Station Locations 
# 1 ( 0, 500, 0) 

# 2 (“433, -250, 0) 

itt (+433, -250, 0) 

//4 ( 0 , 0 , 0 ) 

Theoretical 

Values 

Values From 
Computer 
Simulation 

o 

X 

a 

y 

a 

z 

a 

X 

a 

y 

O 

z 

1. Three station (trilatera- 
tion) solution (# 1 » # 2 , # 4 ) 

4.21 

1.28 

9.9 

4.33 

1.77 

9.9 

2 s Li.n 0 Hrl. 20 d minimum 
variance solution 

D 


Q Q 

y • 

1 Q9 

X 1 ^ 4. 

1.77 

in ^ 

4.V • -V 

3. Explicit minimum variance 
solution 

1.98 

1.31 

6.2 

2.12 

1.65 

7.95 

4. Iterative solution 

1.84 

1.25 

6.0 

1.85 

1.94 

6.61 

1 
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From the above comparisons it may be concluded that the iterative 
solution always provides the lowest value of coordinate standard deviations. 
This is as would be expected. However, the iterative solution causes 
calculation difficulties in that several iterations must be made each 
calculation increment. The linearized minimum variance solutions provide 
good values of standard deviation; however, separate solutions are required 
for the case where all ground stations are in the plane z = 0, and the case 
where ground stations are elevated. Standard deviations for the explicit 
minimum variance solution are close to that of the iterative solution, and 
since the solution is completely general and straightforward, the 
advantages may outweigh the slight disadvantage of increased error standard 
deviations. 


4.3 Position Error Contours 

Error contours have been generated for particular cases to indicate 
the geometrical regions where equal errors occur. These error contours 
have been generated from the theoretical standard deviations of the minimum 
variance iterative solution. 

Figure 4-2 shows a contour plot of the errors in the x coordinate at 
a constant z plane (z = 100 ft). For this and all plots, it has been 
assumed that the two-way slant range measurements have equal random error 
variances, and the plot gives the ratio of the coordinate standard deviation 
to the two-way range measurement standard deviation. In this particular 
figure, contours showing error multiplication factors of unity and two are 
the only integer contours which appear on the plot within the region plotted 
(-1000 to +1000 ft in x and y) . 

Figure 4-3 shows a similar plot of the x coordinate error, except at an 
elevation of z = 500 ft. Figure 4-4 is similar except for an elevation of 
z = 1000 ft. Figures 4-5, 4-6 and 4-7 show similar plots for the y coordinate 
error. The y coordinate error appears comparable to the x coordinate error 
at all elevation planes. 

Error contours for the z coordinate are shown in Figures 4-8, 4-9 and 
4-10. The z coordinate errors, as indicated in Tables 4-1 through 4-6, are 
much more severe than those in x and y. At higher elevation angles, however, 
as indicated in the z = 500 ft plane, the z axis errors become considerably 
smaller . 
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Fig. 4-3. Plot of normalized equal rms error contours in the x coordinate 
(o x /or) in the plane z = 500 ft. The ground station locations 
are shown on the plot. 
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Fig. 4-8. Plot of normalized equal rms error contours in the z coordinate 
(a z /oR> in the plane z = 100 ft. The ground station locations 
are shown on the plot. 
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Fig. 4-9. Plot of normalized equal rms error contours in the z coordinate 
(°z/°r) the plane z = 500 ft. The ground station locations 
are shown on the plot. 
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Figures 4-11 and 4-12 show the effect on the position error contours 
when two of the ground stations (#2 and #3) are elevated 100 ft. Comparing 
Figure 4-11 with Figure 4-8 we see significant improvement in the z axis 
accuracy. The improvement has a larger effect at positions distant from 
the ground stations. Comparison of Figures 4-12 and 4-9 leads to the same 
conclusion. In Figure 4-12, the area of the contours has increased from 
those in Figure 4-8, with the outer contours increasing more than the 
inner contours. 

Figures 4-13 and 4-14 show the effect of changing the ground station 
location to a configuration that favors approaches from the positive y axis. 
The ground station location for this new configuration is shown on the plot. 
All stations are in the plane z = 0 for this configuration. Notice that 
the regions of good accuracy are moved out along the negative y axis, 
whereas the accuracy near the landing point (0,0) is not greatly affected. 

These error contour plots showing the effects of elevating stations 
and moving the ground location of the stations indicate the flexibility 
of the system in tailoring the accuracy to particular approach paths. 

The situation is somewhat analogous to that of obtaining directivity in 
the proper direction with a series of antenna dipoles. If approaches 
will be from a preferred direction, the accuracy in that direction can be 
improved by tailoring the ground station locations accordingly. 

4.4 Position and Rate Errors for Specific Aircraft Trajectories 

4.4.1 General discussion .— As previously discussed, an analysis of 
the coordinate rate errors requires definition of specific trajectories 
because the coordinate rate errors are dependent on aircraft velocity and 
heading, as well as aircraft position. In order to investigate the 
coordinate rate errors, specific straight-in approach trajectories have 
been defined with glide slope angles of 1.3 and 15 degrees. The velocity 
and acceleration profiles of the trajectories are defined in the following 
section. 

The trajectory with the glide slope angle of 1.3 degrees is not within 
the design range of the system, but was included to Investigate the 
severity of errors at extreme low angle approaches. 
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Fig. 4-12. Plot of normalized equal rms error contours in the z coordinate 
(<J z /a R ) in the plane z - 500 ft. The ground station locations 
are shown on the plot {stations #2 and #3 elevated 100 ft). 
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Fig. 4-14. Plot of normalized equal rms error contours in the z coordinate 
(o z /o R ) in the plane z = 500 ft. The ground station locations 
are shown on the plot. 
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4.4.2 Measurement errors In coordinate positions and coordinate rates 
for a straight- in approach trajectory with glide slope angles of 1.3 and 
15 degrees .— This section provides calculations of errors in coordinate 
positions and rates for the four-station multilateration system using a 
straight-in trajectory similar to the Wall St. approach discussed in 
ref. 2. The placement of the ground receivers was varied to determine the 
variation of coordinate errors as a function of ground station location. 

The trajectories considered are shown in Fig. 4-15. The top curve plots the 

velocity profile as a function of distance along the y axis. The velocity is 

constant until a distance of 553 ft from the touchdown point is reached, at 

2 

which point a deceleration at a constant rate of 6.15 ft/sec takes place. 

The final velocity is 1.7 ft/sec at the touchdown point. A plot of altitude 
vs. distance along the y axis is shown in the bottom plot of Fig. 4-15. For 
comparative purposes, two glide slopes are used in the calculations; a 
1.3° glide slope, which represents the Wall St. trajectory, and 15°, which is 
representative of a higher angle approach. The final altitude at the touch- 
down point is 10 ft in both cases. Note that the low-angle glide slope only 
reaches an altitude of approximately 50 ft at 1800 ft from the touchdown point. 

Figure 4-16 provides a sketch of the receiver locations for four different 
cases. The first case la a symmetrical configuration with the receivers 
located on a 500-ft radius circle from the origin of coordinates. Case 2 is 
an L-shaped configuration with three stations located along the axis of the 
trajectory and one station offset along the y axis for a distance of 500 ft. 

Case 3 is a staggered situation with the four receivers offset slightly from 
the ground projection of the trajectory. Case 4 is a diamond-shaped configura- 
tion with one ground receiver placed 1000 ft out along the y axis. In the 
sketches, the location of the transmitter is located by the circled dots. 

In all cases, the ground projection of the trajectory lies along the y axis 
and the touchdown point is at the origin of the coordinates (station 4 
location) . 

The error calculations are made in accordance with the techniques outlined 
in Section 3.4. In all cases, a standard deviation in range measurement 
of 2 ft and a standard deviation in range rate measurement of .1 ft/sec are 
assumed. Plots are then made of the standard deviation of the errors in the 
x, y, and z coordinates, and the corresponding standard deviation of errors 
in the x, y, z coordinate rates. Figures 4-17 through 4-20 show the error 
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Fig. 4-15. Velocity and altitude plotted vs. distance along the y axis 
for assumed trajectories. 



y 


2 ( 0 , 1000 ) 



CASE 2 


- 3 (0,1000) 

* 


<S> i 

(-400,500) 


• 2 

(400,500) 


x 


CASE 4 


:ations for various cases considered, 
lg the y axis. The transmitter is 


63 




S tandar d Deviation , f t 


(> 73 ) 

70 L 


6£> I- 


RG S.O. 2.0 VEL S O. 6 1 


58 


40 


30 


. Case: 1 

Station Coordinates: 

#1 (0,500) 

#2 (-433,-250) 

•a #3 (433,-250) 

#4 ( 0 , 0 ) 

Glide Slope Angle: 1.3° 


20 k 


10 



400 60 © 880 1090 1200 1400 1680 1880 

Distance Along Y Axis, ft 


Fig. 4-17. Standard deviation of errors in coordinate positions plotted vs. 

ground projection of slant range for the specified trajectory. 
Dots are .5 sec apart. Measurement errors are 2 ft rms in 
range and .1 fps rms in range rate. 
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plots for the four cases considered for a glide slope angle of 1.3°. As may 
be seen, the standard deviation in the altitude coordinate (z) is larger than 
that for x and y, as would be expected. 

Figures 4-21 through 4-24 plot the standard deviation in the coordinate rate 
errors as a function of distance along the y axis for the various cases 
considered. The standard deviation of z has a rather large variation in 
most cases, with several peaks and valleys. In all cases, however, the 
standard deviations become small at the touchdown point. 

Figures 4-25 through 4-32 plot similar curves except for a glide slope angle 
of 15°. This glide slope angle is more consistent with the design objectives 
of the multilateration system. For all station locations considered, the 
coordinate position errors remain small (approximately 5 ft rms) , and the 
coordinate rate errors are under 1 ft/sec for the most part. For this 
glide slope angle, the errors are not a strong function of receiver location; 
however, it is possible to attain better accuracy over certain portions of 
the trajectory by judicious location of the ground receivers. 

4.4.3 Measurement errors in coordinate position and coordinate rates 
using elevated ground stations .— Figures 4-33 through 4-40 plot the standard 
deviation of errors in range and range rate using receiver locations 
similar to those shown in Fig. 4-16, except that certain stations are 
elevated 100 ft. The receiver locations for the various cases are 
specified on the figures. For all cases, a 1.3° glide slope is used and 
the basic measurement errors are 2 ft rms in range and .1 ft /sec rms in 
range rate. 

Elevation of one or more of the receivers improves the position 
measurement accuracy significantly in the z coordinate. The major improve- 
ment is achieved at ranges outside the region of station location (greater 
than 1000 ft from touchdown). For example, compare Figure 4-19 with Figure 4-35. 
Measurements of coordinate rates are also more accurate at longer ranges 
with elevated receivers. However, there are locations along the trajectory 
at which large range rate errors occur. For example, see Figure 4-39, where 
the errors peak at a y axis distance of approximately 650 ft. 
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Fig, 4-18. Standard deviation of errors in coordinate positions plotted vs. 

ground projection of slant range for the specified trajectory. 
Dots are .5 sec apart. Measurement errors are 2 ft rms in 
range and .1 fps rms in range rate. 
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Fig. 4-19. Standard deviation of errors in coordinate positions plotted vs. 

ground projection of slant range for the specified trajectory. 
Dots are *5 sec apart* Measurement errors are 2 ft rms in 
range and .1 fps rms in range rate. 
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Fig. 4-20. Standard deviation of errors in coordinate positions plotted vs. 

ground projection ol slant range for the specified trajectory. 
Dots are .5 sec apart. Measurement errors are 2 ft rms in 
range and .1 fps rms in range rate. 
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Fig. 4-22. . Standard deviation of errors in coordinate rates plotted vs. 

ground projection of slant range for the specified trajectory. 
Dots are .5 sec apart. Measurement errors are 2 ft rms in 
range and .1 fps rms in range rate. 
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Fig. 4-23. Standard deviation of errors in coordinate rates plotted vs. 

ground projection of slant range for the specified trajectory. 
Dots are .5 sec apart. Measurement errors are 2 ft rms in 
range and .1 fps rms in range rate. 


71 



S tandard Deviation , f t / sec 


RG S .D 


2 0 UEL S O 


6 


a. 

z 


0 1 


Case: 4 

Station Coordinates: 

#1 (-400,500) 
n (400,500) 

# 3 ( 0 , 1000 ) 

tf 4 (0,0) 

Glide Slope Angle: 1.3 e 


o. 

X 




UJ 


• • 


-L. 


-J- 


O. 


a. 

-I 2L 


>rin 400 600 800 1000 1200 1400 

Distance Along Y Axis, ft 


1600 


1800 


Fig. 4-24. Standard deviation of errors in coordinate rates plotted vs. 

ground projection of slant range for the specified trajectory. 
Dots are .5 sec apart. Measurement errors are 2 ft rms in 
range and .1 fps rms in range rate. 
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Fig. 4-25. Standard deviation of errors in coordinate positions plotted vs 
ground projection of slant range for the specified trajectory. 
Dots are .5 sec apart. Measurement errors are 2 ft rms in 
range and .1 fps in range rate. 
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Fig. 4-26. Standard deviation of errors in coordinate positions plotted vs. 

ground projection of slant range for the specified trajectory. 
Dots are .5 sec apart. Measurement errors are 2 ft rms in 
range and .1 fps rms in range rate. 
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Fig. 4-27. Standard deviation of errors in coordinate positions plotted vs. 

ground projection of slant range for the specified trajectory. 
Dots are .5 sec apart. Measurement errors are 2 ft rms in 
range and .1 fps rms in range rate. 
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Fig. 4-28. Standard deviation of errors in coordinate positions plotted vs. 

ground projection of slant range for the specified trajectory. 
Dots are .5 sec apart. Measurement errors are 2 ft rms in 
range and .1 fps rms in range rate. 
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Fig. 4-29. Standard deviation of errors in coordinate rates plotted vs. 

ground projection of slant range for the specified trajectory. 
Dots are .5 sec apart. Measurement errors are 2 ft rms in 
range and .1 fps rms in range rate. 
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Fig. 4*30. Standard deviation of errors in coordinate rates plotted vs. 

ground projection of slant range for the specified trajectory. 
Dots are .5 sec apart. Measurement errors are 2 ft rms in 
range and .1 fps rms in range rate. 
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Fig. 4-31. Standard deviation of errors in coordinate rates plotted vs. 

ground projection of slant range for the specified trajectory. 
Dots are .5 sec apart. Measurement errors are 2 ft rms in 
range and .1 fps rms in range rate. 
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Fig. 4-32. Standard deviation of errors in coordinate rates plotted vs. 

ground projection of slant range for the specified trajectory. 
Dots are .5 sec apart. Measurement errors are 2 ft rms in 
range and .1 fps rms in range rate. 
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Fig. 4-33. Standard deviation of errors in coordinate positions plotted vs. 

ground projection of slant range for the specified trajectory. 
Dots are .5 sec apart. Measurement errors are 2 ft rms in 
range and .1 fps rms in range rate. 
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Fig. 4-34. Standard deviation of errors in coordinate positions plotted vs. 

ground projection of slant range for the specified trajectory. 
Dots are .5 sec apart. Measurement errors are 2 ft rms in 
range and .1 fps rms in range rate. 
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Standard Deviation , f t 



Fig. 4-35. Standard deviation of errors in coordinate positions plotted vs. 

ground projection of slant range for the specified trajectory. 
Dots are .5 sec apart. Measurement errors are 2 ft rms in 
range and .1 fps rms in range rate. 
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Fig, 4-36, Standard deviation of errors in coordinate positions plotted vs. 

ground projection of slant range for the specified trajectory. 
Dots are .5 sec apart. Measurement errors are 2 ft rma in 
range and .1 fps rma in range rate. 
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Fig. 4-37. Standard deviation of errors in coordinate rates plotted vs. 

ground projection of slant range for the specified trajectory. 
Dots are .5 sec apart. Measurement errors are 2 ft rms in 
range and .1 fps rms in range rate. 
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Fig. 4-38. Standard deviation of errors in coordinate rates plotted vs. 

ground projection of slant. range for the specified trajectory. 
Dots are .5 sec apart. Measurement errors are 2 ft rms in 
range and .1 fps rms in range rate. 
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Fig. 4-39. Standard deviation of errors in coordinate rates plotted vs. 

ground projection of slant range for the specified trajectory. 
Dots are .5 sec apart. Measurement errors are 2 ft rms in 
range and .1 fps rms in range rate. 
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Fig. 4-40. Standard deviation of errors in coordinate rates. plotted vs. 

ground projection of slant range for the specified trajectory 
Dots are .5 sec apart. Measurement errors are 2 ft rms in 
range and .1 fps rms in range rate. 





4.4.4 Summary .— It should be kept in mind that the curves presented 
are for a specific trajectory as given in Fig. 4-15, and that the error 
characteristics are a strong function of the trajectory. It may be 
concluded, however, that within the ranges considered (up to 1800 ft), 
and with glide slopes on the order of 15° , the error characteristics are 
not strongly dependent on specific receiver locations on the ground. For 
glide slope angles of 1-2° , the ground receiver location can 
be used to improve accuracy in the z coordinate. The best ground station 
configuration for the flat trajectory of this type appears to be Case 3, 
in which the ground receivers are staggered along the path of the trajectory. 
Using this configuration, it was possible to maintain the standard deviation 
in the z coordinate near 10 ft for the last 700 ft of the trajectory. It is 
possible that other ground station configurations (or more ground receivers) 
will improve the z measurement accuracy for flat trajectories, and this 
will be investigated in future work. 


4.5 Errors in Estimation of Uncompensated Time Delay Bias 

The technique for removal of unknown time delays in the two-way path 
is described in Section 3.7. Numerical calculations have been made of the 
estimation errors to be expected, assuming initial acquisition of the 
aircraft at various points. The results of these calculations are shown 
in Table 4-7, which gives the initial standard deviation of the time 
delay bias estimation error at four coordinate points. Table 4-8 provides 
the same information except that two of the stations (no. 2 and no. 3) have 
been elevated 100 ft. 

As may be seen, the numerical values of the standard deviation of the 
errors in the bias removal change considerably depending on the coordinate 
point at which initial acquisition takes place. However, the values are 
generally under 3 ft rms on a single calculation basis. Using optimal 
low-pass filtering with a minimum equivalent time constant of 10 sec 
will reduce these initial standard deviations to 1/10 of the values 
shown in the table. If a longer equivalent time constant proves feasible, 
additional reduction in the standard deviation may be achieved. Figure 4-41 
is a plot of the reduction in standard deviation of the time delay 
bias estimation error plotted vs. time, assuming a 10 sample/sec data rate. 
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Table 4- 7 # 


Initial standard deviation of time delay bias estimation 
error at various coordinate points (all stations in z = I 
plane) . 



Station Locations: 
(ft) 


1/1 ( 0 , 500 , 0 ) 

//2 (- 433 ,- 250 , 0 ) 

#3 ( 433 ,- 250 , 0 ) 

//4 ( 0 , 0 , 0 ) 
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Table 4_8. Initial standard deviation of time delay bias estimation 

error at various coordinate points (two elevated stations) 


Point (ft) 


(1000,1000,100) 


(1000,1000,500) 


(500,500,100) 


(100,100,50) 


Station Locations: 
(ft) 


Initial Standard Deviation of Error 
in Paths Specified (ft) 

2R- 

R- + R 0 

R. + R 0 

R-. + R/ 

1 

1 2 

1 3 

1 4 

1.08 

.24 

.78 

1.91 

.49 

.12 

.49 

2.90 

.59 

.02 

1.24 

2.16 

.17 

.73 

3.09 

.006 


//I (0,500,0) 

H2 (-433,-250,100) 
in (433,-250,100) 

a 4 (o, o, o) 




























NORMALIZED STANDARD DEVIATION 




The initial standard deviation drops quite fast during the first second and 
then slowly approaches an asymptotic value of .1 in approximately 10 sec. 

From the above calculations, we may conclude that with four or more 
ground stations, it is feasible to determine the two-way uncompensated time 
delay bias to an accuracy of .3 ft using an optimal filter with an equiva- 
lent time constant of 10 sec. This conclusion is based upon an assumed 
measurement error in the two-way path of 2 ft rms. 
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5.0 COMPUTATIONAL REQUIREMENTS FOR THE MULTILATERATION SYSTEM 


5 . 1 General Requirements 

This section describes an investigation into the computer hardware 
and software requirements for the multilateration system. The basic goal 
is to solve for position and rate ten times per second, using data from 
four or more stations. This should be done on a minicomputer system costing 
approximately $10,000. Use of the weighted least squares algorithms 

described in Section 3.0 is assumed. These require floating point 
arithmetic. 

The input to the minicomputer will consist of 16-bit digital data 
from the four ground receivers. It is anticipated that the ground 
receiver data will be multiplexed so that the received data will consist 
of 8 16-bit words for each measurement cycle. These words will consist of 
the four range measurement words and four range rate measurement words 
(or counts) to be transmitted as fast as possible during each measurement 
cycle (.1 sec). Output from the computer system will consist of a 16-bit 
digital line plus a data ready line for providing output to the navigation 
computer via a data link. In addition to the digital 1/0 channel, it is 
necessary to provide analog outputs to drive x-y plotters and strip chart 
recorders. Eight channels should be sufficient for monitoring the 
operation of the system. The eight analog output channels could be 
utilized as follows: 

Output Analog Channels and Function 


1 . 

x coordinate 

output 

2. 

y coordinate 

output 

3. 

z coordinate 

output 

4. 

x coordinate 

output 

5. 

y coordinate 

output 

6. 

z coordinate 

output 

7. 

ground projection of slant range 

8. 

slant range. 



For convenience in programming the minicomputer, a teletype system and 
high-speed paper tape reader are desirable as input devices. Computer 
timing requirements are considered in the following section. 
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5.2 Timing Tests on a Representative Minicomputer (PDP-8/E) 

To determine rough estimates of computation time, timing measurements 
were made on the weighted least squares computation of position only, 
using the following equation: 

X = (P T P)- 1 P T Q (5-1) 

where P is a 4 x 3 constant matrix, Q is a 4 x 1 vector of modified measure- 
ments and 



T \l> 


AA 


Here, 


" A x 0 0 0 

0 A 2 0 0 

0 0 A 3 0 

0 0 n a 

* 4 


and 


rl> 


AA 



1 

-1 

-1 

-1 


-1 -1 - l " 

5 11 

15 1 

115 


(5-2) 


(5-3) 


(5-4) 


Some initial simplifications in this algorithm were made to reduce the 
required computation. Note that 



(5-5) 
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where 



was also found to have a relatively simple closed form: 



(5-6) 


(5-7) 



Note also that the variance cancels out in (5-1) , so it need not be 
included in the calculation. 

A Fortran program was written to solve for x with the following steps 
repeated a specified number of times in a DO loop to allow timing with the 
sweep second hand of a clock. 
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1) Using precomputed values of through A^, compute the non-zero 


entries in the matrix, equation (5-8) . 


T -1 


2) Call the matrix multiplication subroutine GMPRlS* to form P . 

3) Call GMPRD*to form P T P. 

4) Call GMPRD*to form P 4 >~q Q. (Computation of the q^'s from the 

A. ’s was not included in the timing test.) 

1 T -1 

5) Call the matrix inversion routine MINV to invert P 4>^q P. 

6) Call GMPRD*to compute x. 

7) Compute z from the equation 


- -2 - 2 . 1/2 
z = (2y - x - y ) 


(5-9) 


When this program was run under the Fortran II system on the PDP-8/E, 
the time per position computation was found to be approximately 0.9 seconds, 
clearly an unacceptable figure and somewhat greater than expected. An 
investigation of the Fortran II software revealed that it was not making use 
of the Extended Arithmetic Element (EAE) , and thus required approximately 900 ys 
to perform a floating point multiplication. A call to the manufacturer revealed 
that no version of Fortran II is available which does employ the EAE. 

Next, the program was run under the new Fortran IV system, which does 
employ the EAE. Quite surprisingly, the time per position computation was 
still almost 0.9 seconds. This was traced to the fact that all variables, 
including integers , are represented and operated upon as floating point 
numbers in PDP-8/E Fortran IV. The gain in speed for each floating point 
operation was offset by the many additional floating point operations for 
array indexing, etc. 

To overcome these difficulties, a third version of this program was written 
in PAL-8 assembly language. The required subroutines GMPRD and MINV (matrix 
inverse) were also manually translated to assembly language. Floating point 
operations were performed by calls to the EAE Floating Point Package sub- 
routines described in Chapter 8 of ref. 3. Time per position calculation 
with this software was improved by a factor of ten to approximately 90 ms. 

This is still not acceptable, since both position and rate must be computed 
in 100 ms. However, this time is within a factor of two of being acceptable, 
so that use of the PDP-8/E might be possible with further algorithm simplifi- 
cation (see next section) . (During the course of using the EAE Floating 
Point Package, which is distinct from the floating ooint routines used by 
the Fortran systems, a "bug" was discovered: If a number less than one in 
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*Matrix product subroutine. 



magnitude and a zero are added or subtracted, unnecessary roundoff error 
results. A patch to the Floating Point Package was written to correct 
this, and a letter sent to the manufacturer informing him of this problem 
and the suggested solution.) 


5.3 Increasing Algorithm Efficiency 

Reference 3 quotes the following "typical times" for floating point 
operations on the PDP-8/E with the EAE Floating Point Package: 


Add 

Subtract 

Multiply 

Divide 


160 ys 
180 ys 
200 ys 

160 ys or 190 ys • 


Since these are almost two orders of magnitude greater than the average 
single instruction time, it is reasonable to compare algorithms 
on the basis of the number of floating point ooerations they entail. 


A count of floating point operations for the position computation 
algorithm described in the previous section is given below: 

Multiplications Additions and 
and Divisions Subtractions 


Direct computation of 

*aq from V s 

12 

0 

P T ip'i via GMPRD 
AQ 

48 

48 

( pT i|T^) P via GMPRD 

36 

36 

( pT ip'*) Q via GMPRD 

12 

12 

T -1 

Inversion of P P 

AQ 

30 

26 

(P T l fT 1 P)" 1 P T Q 

via GMPRD 

9 

9 

A 

Computation of z from eq. (5-9) 
(Square root counted as 7 
multiplications based on 
times given in ref. 2) 

10 

2 


TOTALS 


157 


133 



If one assumes an average time of 180 ys per floating point operation, 

then the total of the 290 operations in this algorithm takes approximately 52 ms. 

The remainder of the 90 ms is taken up by various other operations. Note 

that the inversion of the 3x3 matrix takes only about 20 % of the total 

time, with matrix multiplications accounting for the bulk of the time. 

The efficiency of this algorithm could be improved in several ways. 

Note first of all that the matrix equation (5-8), has several zero entries, 

and that the matrix P (see Section 3.0) has one column of 1 ' s . Taking advantage 

of these facts reduces the number of multiplications and additions required 

T -1 

in forming the various matrix products. Also the matrix P P is symmetric, 

so that only six of its elements, rather than all nine, need be computed. 

Finally, it is more efficient to solve a system of linear equations by a 
method such as Gaussian elimination than by inverting the matrix and then 
multiplying by the right hand side vector. The Gaussian elimination algorithm 
is also simpler when the matrix is symmetric, as in the case here. 

An algorithm employing the above ideas is described below and the number 
of required floating point operations tabulated. It involves direct computa- 
tion of the elements of the matrix 


2 r, T ,-1 

°A P % 


-HMl 

1L A 1 A 2 A 3 A 4 J A 2 [_ A 1 2 J A 3 L 1 3J 4l_l 4J 

1_ + + I3 + M _l.ril.l3l 

A 1 |_ A 1 A 2 A 3 A 4 J A 2 L A 1 A 2 J A 3 L A 1 A 3j A 4 |_ A 1 A 4J 

A 1 [ A 1 A 2 A 3 A 4 J A 2 L A 1 A 2 J A 3 j_ A l A 3 J A 4 |_ A 1 A 4_ 


h L h 

a i L A i 


y 2 y 3 

+—+—+— 
A 2 A 3 A 4 


As before, it is assumed that both the an ^ the A^'s are available. 


100 



Multiplications Additions and 

and Divisions Subtractions 


Compute — 
A i 



i = 1, 4 


12 


0 


Compute directly each element 

of °A pT ^AQ * e q uation (5-10) 15 

Compute o 2 ^ P T Q 

via conventional matrix 12 

multiplication 


18 


9 


Compute the six distinct elements 
of the symmetric matrix 
2 T -1 

a A P i P, taking advantage 
of the column of l's in P to reduce 

the number of multiplications 12 18 

Solve for x via symmetrical 

Gaussian elimination algorithm 17 12 

Compute z from (5-9) 

(Square root counted as 7 

multiplications) 10 2 


TOTALS 


78 


59 


Note that the number of floating point operations required, 137, is 
approximately one-half of that in the algorithm used for timing tests. A 
quick look at the weighted least squares algorithm for calculation of rates 
indicated that it would require approximately the same number of operations 
as the above position algorithm. It appears, then, that it is probably 
possible to meet the required 100 ms. time limit for both position and rate 
calculations, with careful assembly language programming. However, there 
would be very little margin for any future program changes or additions. 
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Some additional time could be saved by not updating the weighting matrix 
for the position calculation at each measurement, but rather, say, at every 
other measurement. However, it does not appear that this is possible for 
the rate calculations. The use of an iterative algorithm is another 
possibility. 


5.4 Conclusions and Recommendations 

It appears that weighted least squares calculation of position and rate 
ten times per second is just marginally possible on the PDP-8/E with the 
Extended Arithmetic Element with careful algorithm pruning and assembly 
language programming. However, this would place a limitation on future 
algorithm expansion, and would require a significant amount of programming 
time to make program modifications. 

Barring the discovery of a significantly faster iterative algorithm, it 
appears that the best solution would be to purchase a minicomputer with a 
floating point processor. Although it is believed that this option for the 

PDP-8 costs approximately $7,000, the cost for a floating point processor 
for the Data General Nova series of computers is only $4000. With this 
option, a floating point arithmetic operation can be performed in roughly 
10 ys. This would allow an ample time margin with the existing algorithms, 
and might even permit the operational program to be written in Fortran. A 
complete system (Nova computer with 8K of memory, teletype, floating point 
processor, and necessary interfaces) should not cost significantly more 
than $14,000. 

Floating point operation times for other minicomputers without a floating 
point processor are hard to obtain. However, it is believed that neither the 
Data General Nova 800 with hardware fixed-point multiplication nor the Honeywell 
716 can improve upon the PDP-8/E times . 
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6.0 RECOMMENDATIONS AND CONCLUSIONS 


6.1 Recommended Solution Techniques 

The accuracy studies of the solution techniques indicate that while 
the iterative solution described in Section 3.4 provides the minimum 
variances in the coordinate errors for all cases, it has the disadvantage 
of requiring several iterations within a computational cycle. Possibilities 
of convergence problems also exist if the initial solution is considerably 
in error. For these reasons, it is recommended that a slight amount of 
accuracy be sacrificed in order to use the completely general explicit 
solution described in Section 3.5. Advantages of the explicit solution 
are as follows: 

1. Completely general station locations can be used. The station 
survey positions are entered into the solution as constants, 
without consideration of whether the stations are in a plane or 
not. Elevated stations are handled with ease. 

2. The solution is straightforward and requires no iterations. 

3. The accuracy of the solution differs only slightly from the 
minimum variance obtained from the iterative solution. 

4. Computation times should be reasonable, since a number of 
the matrices and matrix multiplications do not need to be 
recalculated in real time. 

For range rate calculations, the minimum variance solution outlined 
in Section 3.6 is straightforward and is recommended for use. For the 
time-delay bias removal, the technique described in Section 3.7 is 
recommended with an optimal low-pass filter as described in Section 3.7.3. 

This technique provides for calculation of the time-delay in each separate 
two-way path, but does not separate transponder delays from transmitter 
and receiver delays. To date, no way of separating the delay into transponder, 
receiver, and transmitter delay has been determined, except by using a 
separate calibrated transponder at a ground location. 

The solution will have the capability of providing coordinate and 
coordinate rate outputs for three stations in case one station is not 
receiving data. Also, provision can be made for utilization of more than 
four ground stations if desired. 
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A flow chart Indicating the recommended calculation scheme is shown 
in Figure 6-1. 


6.2 Recommended Computer Specifications 

In Section 5.0, the general considerations for the system minicomputer 
were discussed. Based on these considerations and known operational 
requirements, it is recommended that a minicomputer with the following 
characteristics be purchased: 

1. A minicomputer with 16,384 words of memory, a memory cycle time 
of about 1.0 ps, and a preferred word length of at least 16 bits; 

2. A floating point processor to perform floating point arithmetic 
operations in 10 to 20 ps each; 

3. A digital I/O channel providing 16 digital lines plus an interrupt 
line for input, and 16 digital lines plus a data ready line for 
output; 

4. Eight (8) analog outputs for plotter drive; 

5. A teletype interface; 

6. A high speed paper tape reader; and 

7. FORTRAN IV software capability. 

Surveys of various manufacturers were made to determine which mini- 
computer system would provide the above characteristics for a minimum price. 
Based on this survey, it appears that a Data General Nova 11/10 system will 
meet the specifications at a minimum price. The total cost (GSA) of the 
Data General System meeting the above specifications is approximately $14,060. 


6.3 Recommended Future Work 

During this study, geometrical errors were determined for trajectories 
consisting of straight-line segments. In future work, more complex trajec- 
tories should be considered, including spiral descents and curved, decelerating 
approaches. In addition to the study of more general trajectories, additional 
study should be made of the effect of station location on the system 
accuracy. Using existing programs, additional error contours can be 
generated for a large number of hypothetical station locations. 
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Fig, 6^1. General flow chart for recommended solution techniques — 
Operational System, 
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Fig. 6-1. Continued. 
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(SEE SECTION 3.6) 



Fig. 6-1. Continued. 
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Fig . 6-1 . Concluded . 




The benefits to be achieved by using more than four ground stations 
should also be investigated. Since ground stations are anticipated to 
be low cost, and the added computational burden is not great, it is felt 
that additional stations may prove worthwhile in cases where the geometrical 
accuracy requires improvement. 

In the studies discussed in this report, a basic measurement error has 
been assumed. This assumption has been based on analyses accomplished by 
LRC personnel, slanted toward the determination of hardware-associated errors. 
Additional analytical work is required to complete some of the hardware 
error studies. In addition, results from field checks of the equipment 
should be incorporated in the error studies to provide firm numbers for the 
basic two-way path measurement errors . 

Additional consideration should be given to the integration of the 
multilateration system into the VALT navigation system. The solutions 
discussed in this report are all point estimation techniques; that is, 
smoothing of data based on the sequential measurements has not been recom- 
mended. This has been done to prevent time lags in the output data that 
may detrimentally affect the operation of the navigation and guidance 
systems on board the aircraft. It is felt, however, that in' some cases, 
smoothing of the output data may be feasible and that the associated time 
lags could be made compatible with the navigation system. For this reason, 
navigation system personnel should define the transient characteristics 
that are necessary in the basic measurement data for compatibility with the 
navigation and guidance algorithms. 

In future work, extensive simulations of the solution technique 
recommended in Section 5.0 should be conducted to assure that there are 
no points in the coverage volume where the solution fails. In addition, 
considerable work is necessary to determine the most efficient algorithms 
for providing the necessary computations in real time. The associated 
software must be developed and checked out for the minicomputer selected 
for use in the system. 
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appendix a 


An Explicit General Solution to the Trilateration Position Location 
Problem in an Arbitrary Coordinate System 

The general solution outlined in this appendix is similar to a 
solution given in algorithmic form in ref. 4, except that matrix techniques 
are used to permit extension to the case of redundant ground stations. 

For three ground stations, three equations are obtained: 


2 

2 

2 



(A-l) 

R i 

= R 

+ pi - 

2R • 

Pi 

2 

2 

2 




R 2 

= R 

+ p 2 - 

2R • 

P2 

(A-2) 

„ 2 

^2 

2 

2R • 



R 3 

* R 

+ p 3 - 

P3 

(A-3) 

where the nomenclature is defined 

in Fig. 

3-1, 

page 14 . 



The solution is started by subtracting (A-l) from both (A-2) and (A-3) 
to obtain: 


2 2 2 2 
r 2 - V - P 2 + Pi 


= + R • (p x - p 2 ) 


2 2 2 2 
V- R 1 -p 3 +Pl 


= + R • (p x - P 3 ) 


(A- 4) 


(A- 5) 


Now let the left-hand side, which contains the measurements and surveyed 
values, be designated as q^ and q^ so that the equations become 

q x = + R • (P x - P 2 ) (A- 6) 

q 2 = + R • (p 1 - p 3 ) . (A- 7) 

These equations can be written in matrix form as: 

Q = [P] x (A- 8) 
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where 


x 

y 

z 


IP] 


X 1 - X 2 y l ~ y 2 2 1 - Z 2 


X 1 ' X 3 


y l - y 3 


Z 1 ' Z 3 


Now partition the matrix [P] and x aa follows: 


Q “ [A B] J p z 


where 


and 


h*i 

[A] = Xl * 2 7l y *\ 
L*1 - X 3 y l " y 3 J 


[B] 


j Z l " Z 2 
) Z 1 - Z - 


(A-9) 


With this partitioning, we have 


Q = A p + Bz 


and for A nonsingular, it is possible to solve for p as 


p=A^Q-A^Bz 


or 


where 


p » CQ - CBz 


[Cl. 


-1 1 f y l “ y 3 y 2 y l 

-w -i[4-4 4-4 


d = ( x i - x 2 ^ y i " yJ ~ ^ x i " ~ yJ - 


3 /^i 


(A-10) 


(A-ll) 

(A-12) 


From the above equation, x and y are found in terms of z as 
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X = 


k l - d l z 


y = k 2 ■ V 

where d L = jr [(yj^ - y 3 )U 1 ~ z 2 > + (y 2 “ y l^ z l " z 3^ 


- [(x 3 - x 1 )(z 1 - z 2 ) + (x x - x 2 )(z 1 - 


z 3 )) 


k, = 


< y i - y,) + jr ( y > - y i> 


k 2 D ^ x 3 ~ x l^ + D ^ X 1 ~ x 2^ 


Now, writing eq. (A-l) in terms of x, y, z, 


2 2 2 2 2 

= x + y +z + P 2 ” 2x^x - 2y^y - 2z^z 


and substituting equations (A-13) and (A-14) gives the equation 


az + bz + c = 0 


where 


2 2 

a = 1 + dj* + d 2 


b = -2(z x + k 1 d 1 + k 2 d 2 - x x d x - y x d 2 ) 
c = k^ 2 + k 2 2 + Pj 2 - R^ 2 - 2xjk^ - 2y^k 2 


Thus, it is possible to solve for z as 

-b ±Vb^ - 4ac 


2a 


and x and y are given by 


x = k^ - d^z 


y = k 2 ~ d 2 z 


The sign of z in (A-17) is chosen to make z positive. 


(A-13) 

(A-14) 


(A-15) 


(A-16) 


(A-17) 


(A- 18) 
(A- 19) 
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APPENDIX B 


MATRIX FORMULATION OF COORDINATE ERRORS 

The Taylor Series expansion for the errors in coordinates can be written 
as a function of the measurement errors AR as 


where for four stations, 


Ax 


Ax = [D] AR 





AR 1 1 

I Ax 



J- 

AR« 

I Ay 

and 

AR = 

z 

AR 

/ Az 


| 

1 

k 4 


(B-l) 


(B-2) 


and [D] is the matrix of partial derivatives: 


[D] 



9x 

9x 

9x 

3x 


9R^ 

9R 2 

9R 3 

9R 4 

s 


iz_ 

i2_ 

iz_ 


9R^ 

9R 2 

9R 3 

9R, 

4 


9z 

9z 

9z 

9z 


aRj^ 

9R 2 

9R 3 

9R 4 


(B-3) 


To obtain the variance-covariance matrix, multiply Ax by its transpose and 
take expected values. 


Ax Ax T = [D] AR AR T [D] T 


ifrr— - E{Ax Ax ^ } = 
Ax 


[D] i|»_[D] 
AR 


(B-4) 


where is the variance-covariance matrix of the range errors (see Appendix C) 
and is the variance-covariance matrix of the coordinate errors. In cases 

where a coordinate transformation is used of the form 

X = [K]x 
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then the covariance matrix in the new coordinate system is related to the 
covariance matrix in the old system as 



APPENDIX C 


DERIVATION OF THE COVARIANCE MATRIX OF THE RANGE MEASUREMENT 
ERRORS FOR THE PLANNED MULTILATERATION SYSTEM 

I 

The planned system uses two-way range measurements from a single 
transmitter to four or more ground receivers. Thus the measurements 
are of the form, for four ground receivers, 

M 1 " 2R 1 

M 2 = R x + R 2 (C-l) 

M 3 “ R 1 + R 3 

m 4 = R x + R 4 


where the are the measurements assuming the transmitter and receiver 
// 1 are colocated. The one-way ranges are given by 



(C-2) 


If it is assumed that all two-way range measurement errors have equal 
2 

variance o , then 


2 

<AR X AR X > = f- 


2 

<AR X AR j > « - f- 


J - 2, 3, 4 


(C-3) 

(C-4) 


119 




Thus, the covariance matrix for the one-way range errors may be 
written as (for four ground stations) 

1 -1 -1 -1 
-15 11 

-115 1 

-1115 


A similar expression is obtained for the range rate errors if it is 
assumed that all two-way range rate measurements have equal variance. 



(C-7) 
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APPENDIX D 


COVARIANCE MATRIX OF COORDINATE ERRORS FOR A LINEARIZED SOLUTION 
FOR THE SPECIAL CASE OF ALL GROUND STATIONS IN A PLANE 


In this solution as given in Section 3.3.2, the solution is given by 


x 


(P T P)- 1 P T 



(D-l) 


where is the covariance matrix of the range difference measurements. 
The covariance matrix of errors in the vector x is found by taking the 
expected value of Ax Ax to obtain: 


where 


rj) — = (p T ijT 1 p) -1 
Mx 11 V AQ ' 


ill as T ill T 

V AQ ^AR 


(D-2) 


(D-3) 


and the matrix T is defined in Section 3.3.2. 
Since in this case, the vector x is 


x = 


x 

y 

u 


2 

where u = R /2, the variance of the z coordinate errors must be found from 
a Taylor Series expansion of the equation for z, This equation is 


z = [2u - x 2 - y 2 ] 1 ^ 2 = F(x, y, u) 


(D-4) 


Thus we have 


Az = [D] Ax 


(D-5) 
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where D is a matrix of partial derivitives given explicitly by 



(D-6) 


This variance of errors in the z coordinate is therefore 



APPENDIX E 


COVARIANCE MATRIX OF COORDINATE ERRORS FOR A LINEARIZED SOLUTION 
FOR THE SPECIAL CASE OF ALL GROUND STATIONS NOT IN A PLANE 

For this solution, as given in Section 3.3.3, we have for four ground 
stations : 

x = (P T P) _1 P T Q . (E-l) 

The covariance matrix for errors in the vector x is 

= (P T P) _1 P T ^ P(P T P) _1 (E-2) 

where is found as shown in Appendix D and 

! x 

y 

z 

For more than four ground stations, the weighted least squares 
technique may be used, and the associated covariance matrix is 



T -1 -1 

if; — = (p th x p) ± 

^Ax v ^AQ } 


(E-3) 
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APPENDIX F 


THE OVER-DETERMINED (ITERATIVE) SOLUTION FOR THE DATA OBTAINED 
FROM THE RF MULTILATERATION SYSTEM 


F.l Over-Determined Solution for Four Ground Stations 


The methods for obtaining the following are given below: 

(1) The over-determined solution for x, y, z — as 
obtained from four direct (basic) measurements, 

A, B, C, and D. 

(2) The variance-covariance matrix for the errors in x, y, z. 
Since the errors in A, B, C, and D are assumed to be uncorrelated 

and to have equal standard deviations (viz 10 ft) , it is believed to be 
advantageous to derive (1) and (2) above from A, B, C, and D directly — 
and hence from the ellipsoids which they represent. 

That is. 


2 >/(x-x 1 ) 2 + (y-yj) 2 + (z-z^ 2 

Vvx-x^ 2 + (y-yj) 2 + (z-z ^) 2 +V(x-x 2 ) 2 + (y-y 2 ) 2 + (z-z 2 ) 2 

V(x-x 1 ) 2 + (y-y^ 2 + (z-z^) 2 + V(x-x 3 ) 2 + (y-y 3 > 2 + (z-z^) 2 

V(x-x 3 ) 2 + (y-y^) 2 + (z-z^) 2 + V(x-x^) 2 + (y-y^) 2 + (z-z^) 2 


= A 1 (F-l) 
= A 2 (F-2) 
= A 3 (F-3) 

= A 4 . (F-4) 


Note: If the attitudes of the transmitter and receivers are zero, then z^ 

z 2 = z 3 - z^ = 0. For "exactness", the survey positions should be used. 

Now expand the left-hand sides of the above equations in a Taylor 
Series expansion around an initial estimate, e.g., the presently available 
solution, neglecting all except the first-order terms. The results are 
as follows: 


f 1 (x,y,z)] o + 


34 il 3£ il 3f il 

a7> 4 * + 3yjo 4y + Fiji 


Az 


A 


(F-5) 
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f 2 ( x »y » z ) ] Q + 



9f 

Ax + r— 1 
o 9y 


i ay + SL 


Az = B 


3f| 9fl 3 f ."I 

i (x -y- z)1 o + ST-Jo ix + arjo 4y + Tri 

9f,l 9 f ,”] 3f /] 

(x,y,z)] + r Ax + Ay + 

r ,J * o 9x Jo 9y Jo J 9z Jc 


Az = C 


Az - D . 


(F-6) 


(F-7) 


In the above, 

f 1 (x,y,z)] Q = 2\/(x-x 1 ) 2 + (y-y^) 2 + (z-z.^ 2 ( p -9) 


is evaluated at (x,y,z) = (x q y^ z q ) , which is the initial (starting) 
estimate for x,y,z. Similarly, f 2 (x,y,z)] Q and the others are so defined. 


Now define the following: 


> 

> 

II 

A x ~ f-, (x,y,z)] Q 




AA 2 = 

A 2 - f 2 (x,y,z)] Q 




aa 3 = 

a 3 ' f 3 ( x *y» a >l 0 



(F-10) 

> 

> 

II 

A 4 - f 4 (x,y,z)] o . 

3f l 

af l 

9 f l 

simplify the notation from 

l 

9x . 

to - — and so on. 
Jo 9x 

^ ’ 37" is 


to be evaluated at (x,y,z) = (x o ,y o ,Z Q ) — ° r at new (i m P rove< l) values 
in subsequent iterations. 

The notation can be further simplified to the following: 


[D] 


4x3 



(F-ll) 


3x1 


4 '4x1 


where 
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Matrix 

of 

Partial 

Differentials 


Hi 

af i 

Hi 

9x 

3y 

9 2 

3£ 2 


af 2 

9x 

8y 

3z 

Hi 

3f, 


9x 

3y 

9 z 


Ha 

af 4 

9x 

dy 

3 z 


For further simplification of notation, let 

A? = A y I 


\ 4/4x1 

05 '“UxJ ' 


Consequently, 


^4x3 ^3x1 ^4x1 


Therefore, 


T T 

D D Ax = D AA , 


D = the transpose of D. 

The least-squares solution for Ax is thus obtained from the following: 


Ax = (D T D)" 1 D T AA , 


-1 T 

- 6 D AA , 
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where 


G = D T D (F-18) 

and 

G ^ = the inverse of G. 

The value so obtained for Ax provides an improvement over the original 
estimate of (x,y,z) = (x o ,y o ,z Q ) anc * t ^ e a ^ ove P rocess is repeated until no 
further change (improvment) in Ax results. 


F.2 Variance-Covariance Matrix of Errors in x, y, z 

In the above expression for Ax, let us multiply both sides by their 
respective transposes and then take the expected values. The result is 
the following: 

V iK * G ' 1 1,1 V 4A D(G ' 1)T ’ 


-1 T -1 

= G D V D G (since G is symmetric), 

tAri. 


(F-19) 


where 


and 


V. the variance-covariance matrix of Ax 
Ax ~ 


V 


AA = the variance-covariance matrix of AA 


= a 


AA 


1 0 
0 1 
0 0 
0 0 


= a 2 I 
AA 


0 

0 

1 

0 


0 

0 

0 

1 


(F-20) 


where 

Hence, 


I = the identity matrix. 


2 -1 T -1 

V - of. G D I D G 
Ax AA 


2 -1 T -1 

o. A G D D • G 
AA 
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(F-21) 



Since 





Ax 




F.3 Calculation of the Partials 


R t = Jcx-xJZ + (y-y^) 2 + (z-z ± ) 2 


for i = 1, 2, 3, 4. 

Consequently, 



(F-22) 


(F-23) 


(F-24) 
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F.4 Advantages of the Over-Determined Solution 

The following are advantages of using the over- determined solution: 

(1) All of the observational data get used in a valid least-squares 
solution. 

(2) It provides the standard (i.e., the best attainable solution) 
against which all simplifications and shortcuts can be compared. 

(3) It indicates whether a particular computational simplification 
is warranted, or whether a more nearly exact solution is essential. 

(4) The variances and also the covariances of the errors in x, y, and 
z are readily obtained via matrix algebra in the calculation of the GDOPs 
associated with the over-determined solution. It is believed that the 
correlation structure (and hence the variance-covariance "picture") may 
change significantly throughout the track of the helicopter. This implies 
that the "uncertainty ellipsoid" associated with x, y, z is changing its 
orientation throughout the track. In fact, if the variance-covariance 
matrix at one point cannot be obtained at another point via the multiplication 
by a single scalar, a case for Kalman filtering can be made. Depending on 
the speed with which the correlation structure changes, Kalman filtering may 
produce spectacular improvements in accuracy. 

(5) The matrix algebra used in obtaining GDOPs for the over-determined 
solution can also readily be used to obtain GDOPs for the exactly determined 
solution. 


F.5 Disadvantages of the Over-Determined Solution 

(1) The major disadvantage is the necessity of using iterative 
techniques with associated loss of computational time. 

(2) Problems of convergence of the solution may arise if the initial 
solution is considerably in error. 
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APPENDIX G 


DERIVATION OF ERROR VARIANCES FOR THE EXPLICIT MINIMUM 
VARIANCE SOLUTION GIVEN IN SECTION 3.5 


In this solution, the variance of the z estimate is given by 
2 T -1 „.-l 

°a£ = (U *Az U) (G-l) 

where U is a n x 1 unit matrix and tfi. is the covariance matrix of the 

r Az 

n z(i) solutions obtained from the equations 


a z (i) + b^ z(i) + c^ = 0 


(G-2) 


To determine the first step is to obtain an expression for the 

variation in z(i) in terms of variations of the coefficients a, b, and c. 
Taking differentials in eq. G-2 and solving for Az(i) gives, since Aa = 0, 

Mi) = (lit l bj 4b i - ^2az 1 + 4c i i - 1.0 • (G - 3 > 

In the terms in parentheses, z is used instead of z(i) as an approximation. 
In matrix form 


where 


Az = [K] [zAb + Ac] 


2az + b. 


[K] = - 


2az + b. 


ETC 


_ln x n 


(G-4) 
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Ab . 


ac ii 


Az(l) v 


Ab, 


AC 2 


Az(2)l 

Ab = i 

• 

• 

Ab 

n 

Ac = 

L x i ^ 

• 

• 

Ac 

n 

, and Az = 

^ n x 1 

• I 

i Az(n) ) 


The differentials Ab^ and Ac^ are found from the expressions for 
and (see Section 3.5) 


b i = - 2z i - 2B T C T CQ + 2p i T CB 
T T 

Ab ± = - 2B C CAQ 


Ab = - 2UB T C T CAQ 


(U = n x 1 unit matrix) 


Ab - - 2UB T C T CTAR 


2 T T T T 2 

C ± = - R^ + P i P i - 2p i i CQ + QVCQ + z t 

AC ± = - 2R i AR ± - 2p t T CAQ + 2Q T C T CAQ 

AC = - 2R AR - 2yC AQ + 2UQ T C T AQ 


AC 

where 


,T„T 


2 {-R + (UQ C C - yC) T} AR 


[R] = 


* R, 

0 


ft 



’ R ! 

O 

1 

& 

1 

1 






1 

n 

0 

• 

R 2 

• 

ft 

ft 

3 


[T] = 

0 

ft 

R 2 ' - R n 

. R , -R 
n-1 n 


• 

ft 

R 



ft 

-R 




n. 

n x n 



n J 


(n x 


Y = 


n 
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J n x 2 


(G-5) 

(G-6) 

(G-7) 

(G-8) 

(G-9) 

(G-10) 

(G-ll) 

(G-12) 


1) x n 



Now, letting 


V = -2UB T C T C (n x n) (G-13) 

and 

S - 2 {-R + (UQ T C T C - yC) T} (n x n) (G- 14) 

we have 

Az - K[z V + S] AR (n x 1) (G-15) 

and 

i|T— = E {Az Az 1 } = K[zV + S] i (r— [zV + S] T K T (G-16) 


where is the covariance matrix of the one-way range measurement errors. 
The covariance matrix of p is found as follows. Since 

p = CQ - CBz 

then 

Ap *» C AQ - CB Az 
Ap « CT AR - CB Az 

but 

Az = (U^ \jr~ U) ^ U T t(/~ Az 
Az = F Az 
therefore 

Ap = CT AR - CBFK [zV + S] AR 
Ap - {CT - CBFK [zV + S]} AR . 

Let 

W - CT - CBFK [zV + S] 

then 

tp - = W t|r— W T , (G-22) 

Ap AR 

where, as previously indicated is the covariance matrix of the one-way 
range measurement errors (see Appendix C) . 


(G-17) 

(G-18) 

(G-19) 

(G-20) 

(G-21) 
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Error Variances when Matrix \ p— cannot be Inverted 
Az 

In some cases, it may be impossible to invert the matrix tfj— - because of 
high correlation between errors in calculation of the z(i). (i.e., roundoff 

errors in computation may preclude the inversion.) In this case, a procedure 
as follows is suggested: 

1. Use the average (unweighted) value of the z(i) to calculate — . 

Az 

2. Select the z estimate z as the value of z(i) which has the minimum 

variance. Variances are determined from the diagonal of i |r— . 

„ ~ Az 

A A A 

3. Calculate x and y as before, using the estimate z in the equation 

for p . 

When the above procedure is followed, the variance of errors in z will be 
as determined from the selected diagonal term of if^-, and the variance of the 
x and y estimate is calculated as follows: 

Since 


Ap = CAQ - CBAz 


(G-23) 


and now, 


AZ 


2a < z> + b 


( < z > 


Ab^ + 


AC i ) 


(G-24) 


where the nomenclature has been previously defined. The variations in 
and are given by (i is index for selected z(i) value) 

Ab ± = - 2B T C T CT AR (G-25) 

AC. = [ - 2 R' - 2p 1 T CT + 2Q T C T CT] AR , (G-26) 

where E* = (a^, a 2> a^, a^} 
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■j- 1 * i - 1 

Sj - o, } + i . 


For simplicity in notation, let 


k 


i 


-1 

2a <z > + 


(G- 27) 


S' - [ - 2R' - 2?^^ + 2Q T C T CT] 


(G-28) 


and 


V - - 2BC T C T 


(G-29) 


Then we have, 

Az = k [ < z > V + S ' ] AR 

and 

Ap = {CT - CBk^ [ < z> V + S']} AR 
or 

Ap = W' AR , 


(G-30) 

(G-31) 

(G-32) 


so that we have for the covariance matrix of errors in x and y estimates, 



tlr-rr W 

y AR 


,T 


(G-33) 
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APPENDIX H 


ERROR COVARIANCE MATRIX FOR COORDINATE RATES 


In Section 
range and range 


3.6, the solution for coordinate rates in terms of the 
rate measurements is given as 


x 


T -1 -IT 

(S * aA s) s 



(H-l) 


where the A are the measurements and the remaining matrices are defined 
in the text. For simplicity, define 


[H] = (S T S)' 1 



(H-2) 


so that 


x « H A 


(H-3) 


Now, recalling that [H] is a function of the position (x, y, z) estimate, 
small error in x is given by 


Ax = HAA + AH A 


(H-4) 


Since A = Sx, eq. (H-4) may be written as 
Ax = H AA + AHSx . 

Also, postmultiplying eq. (H-2) by S gives 
HS = I 

where I is an identity matrix. Thus 
AHS + HAS = 0 


(H-5) 


(H-6) 


(H-7) 
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or 


AHS = - HAS 


(H-8) 


Therefore, 


Ax = H AA - HAS x 


A JL 

and since x « Ax/At we can also write 


(H-9) 


x — — > _ 

Ax = H AA - HS Ax 


(H-10) 


The covariance matrix of the errors is found by taking the expected 
— — x 

value of Ax Ax as follows : 


*Ax = E A ^ T } " 


H 


*AA H 


+ HS 




’T T 
S H 


(H-ll) 


— -X 

where it is assumed that the covariance terms E{AA Ax } are zero (i.e., no 
correlation between range rate measurements and position estimates). 

The first term in eq. (H-ll) can be reduced to obtain for the coordinate 
rate error covariance matrix, 



T — 1 — 1 • • T X 

(S $ • S) + HS ill - S H 
AA Ax 


(H-12) 


where is the covariance matrix of the coordinate position errors and the 
remaining matrices are defined above. 
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